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Abstract 

The northwest Australian summer monsoon owes a notable degree of its interannual 
variability to interactions with other regional monsoon systems. Therefore, changes in 
the nature of these relationships may contribute to variability in monsoon strength 
over longer time scales. Previous attempts to evaluate how proxy records from the 
Indonesian-Australian monsoon region correspond to other records from the Indian and 
East Asian monsoon regions, as well as to El Nino-related proxy records, have been 
qualitiative, relying on ‘curve-fitting’ methods. Here, we seek a quantitative approach 
for identifying coupling relationships between paleoclimate proxy records, employing 
statistical techniques to compute the interdependence of two paleoclimate time se¬ 
ries. We verify the use of complex networks to identify coupling relationships between 
modern climate indices. This method is then extended to a set of paleoclimate proxy 
records from the Asian, Australasian and South American regions spanning the past 
9,000 years. The resulting networks demonstrate the existence of coupling relation¬ 
ships between regional monsoon systems on millennial time scales, but also highlight 
the transient nature of teleconnections during this period. In the context of the north¬ 
west Australian summer monsoon, we recognise a shift in coupling relationships from 
strong interhemispheric links with East Asian and ITCZ-related proxy records in the 
mid-Holocene to significantly weaker coupling in the later Holocene. Although the 
identified links cannot explain the underlying physical processes leading to coupling 
between regional monsoon systems, this method provides a step towards understand¬ 
ing the role that changes in teleconnections play in millennial- to orbital-scale climate 
variability. 

“Corresponding author: fiona.mcrobie@research.uwa.edu.au 
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1 Introduction 


The northwest Australian summer monsoon, and the related circulation over the Mar¬ 
itime Continent (i.e. the Indonesian-Australian summer monsoon - IASM), is a critical 
feature of the global low latitude circulation. It provides a global heat source, and is the 
primary region of latent heat release associated with both the Southern Oscillation and 
the Madden-Julien Oscillation (MJO; McBride 1998; Hung and Yanai, 2004). Despite its 
importance, the Australian summer monsoon, occurring over the northwest Kimberley re¬ 
gion of Australia, is relatively shallow, with sensible heating only observed below 750 hPa 
(Hung and Yanai, 2004). Monsoon precipitation is relatively low, with annual November 
to April precipitation over northwestern Australia ranging from a mean of 1200mm (Kim- 


berley Coastal Camp; 

Bureau of Meteorology 

2014b 

) in the northwest, to 500nmi at the 

south (Jubilee Downs, 

Broome; Bureau of Meteorology 

2014b 

), over a distance of some 


500km. Such a relatively weak monsoon system, located at the southern margins of the 
more general IASM regime, should be sensitive to changes in forcing mechanisms acting at 
both the global and regional scale, and over short and long time scales. 


While a range of considerations come into play (e.g. Chang et al., 1979; Hung and Yanai 


2004; Wheeler et al. 2009), the dominant control on the Australian summer monsoon 


relates to the controlling role of the thermal land-sea contrast that manifests itself in 
the heat lows that develop during the summer months. IASM strength is also tied to 
the latitudinal position of the Intertropical Convergence Zone (ITCZ), separating equator- 
ward easterlies from poleward westerlies. The monsoon regime is characterised by summer 
rainfall associated with low-level westerlies that extend from the equator to around 15°S. 
The position of these westerlies is associated with the monsoon trough, representing a 
broad zone of strong convective activity with generally westerly inflow and characterised 
by the occurrence of monsoon depressions and tropical cyclones, defining the southern edge 
of the IASM region. With the progression of the seasons there is a northward displacement 
of the ITCZ, such that by the boreal summer it is located well to the north of the Maritime 


Continent, and is now associated with the East Asian summer monsoon (Chen et al., 2004) 


It is the onset of westerly flow which defines the Australian summer monsoon circula¬ 
tion, and ‘active’ monsoon phases are linked to the MJO, resulting in strong convective 


activity and precipitation over the monsoon region (Hung and Yanai 2004 Wheeler et al. 


2009). Interhemispheric interactions between the IASM and the Northern Hemisphere are 


provided by cold surges emanating directly out of the East Asian winter monsoon, and 
leading to strong convective activity in the South China Sea and over the wider IASM 
region (Chang et al., 1979). It has also been suggested that these cold surges may also 
be directed into the Arabian Sea, enhancing MJO activity (Wang et al. 2012a), which 


provides a link with the Northern Hemisphere. These relationships make it clear that the 
present IASM is driven by an ensemble of regional and global scale climate controls (e.g. 


Chang et al. 

1979; 

Meehl 

1987 

Hung et al. 

2004 

Wang et al. 

2012a) 


When considered over longer time scales, additional drivers at both the global and 
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regional scale need to be introduced. Milankovich insolation forcing of global monsoon sys¬ 
tems has been long recognised (e.g. Clemens et al.| 1991 Bowler et al., 2001 Wang et al. 
2008). Coupled ocean-atmospheric modelling studies have sought to explain the response 


of the northwest Australian monsoon to direct insolation forcing (Liu et al., 2003 Wyrwoll 


et al. 2007, 2012). These results suggest that although precession dominates changes in 


Northern Hemisphere monsoon strength, the Australian monsoon response is also signifi¬ 
cantly impacted by ocean temperature feedbacks (|Liu et ah 2003) and tilt forcing) Wyrwoll 


et al., 2007). Liu et al. (2003) suggest that the enhanced Australian monsoon at 11,000 


years BP, contrary to reduced summer insolation, is due to a combination of sea surface 
temperature feedbacks and inflows from a strong East Asian winter monsoon. 

The interconnected nature of these coupling relationships provides evidence for the 


‘global monsoon’ model as advocated in recent literature (Trenberth et al. 2000; Wang 


et al., 

2009, 

2012b, 

2014 


This concept has been advanced to portray monsoon activity 
as a single body of tropical convection migrating about the equator according to seasonal 
heating, and tied closely to the positioning of the ITCZ (Wang et al. 2009, 2014). Over 


longer time scales, a coherent response of regional monsoons to Milankovich insolation 


forcing is noted by Kutzbach et al. (2008). Using an accelerated transient simulation span¬ 


ning 284,000 years, the authors display a positive response in regional monsoon systems to 
orbital forcing, with lead/lag relationships driven by local land and sea surface tempera¬ 
ture feedbacks. As such, the global monsoon model has been extended to the paleoclimate 


context to describe this somewhat synchronous response to orbital forcing (Ziegler et al. 


as well as abrupt events such as the Heinrich Stadials (Cheng et al. 2012). 


Here, we use complex network theory to analyse relationships between the northwest 
Australian summer monsoon, related monsoon systems and likely forcing climate states. 
We explore these relationships within the context of the ‘global monsoon’, and through 
this we seek to separate global, interconnected relationships and drivers from more local 
controls. Using this approach, we attempt to establish the changing nature of the dynamical 
coupling relationships of the Australian summer monsoon over Holocene time scales. 


2 Methods 

Complex network theory offers a method for identifying coupling relationships and long- 
range teleconnections by connecting ‘similar’ data sets. As such, it provides a suitable 
approach to assess interactions between monsoon systems within the context of the global 
monsoon (Donges et al. 2009). By defining a measure of similarity between climate 


time series, climate networks have been shown to provide insight into dynamical inter¬ 
actions beyond the scope of traditional statistical analysis (e.g. Donges et al. 2009 


2013; van der Mheen et al. 2013, Peron et al., 2014). Measures of similarity include 


linear cross-correlation, mutual information, and event synchronisation between extremes 
(Donges et al., 2009 Rehfeld and Kurths, 2014). Applying complex network methods 
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to modern climate data is relatively straightforward, due to the availability of gridded 
datasets and high-density observation networks, but they also provide a powerful tech¬ 
nique for analysing paleoclimate time series. This is demonstrated by Rehfeld et al. (2013) 


who developed a paleoclimate network of the Indian and East Asian summer monsoons 
covering the past 1,100 years, demonstrating distinct changes in network structure be¬ 
tween the Medieval Warm Period, Little Ice Age and present day. The application of these 


techniques is facilitated by the development of a Matlab toolbox (Rehfeld and Kurths 
2014; http://tocsy.pik-potsdam.de/nest.php). Here, we first construct a climate net¬ 


work using modern convective indices to demonstrate the veracity of complex network 
theory to identify dynamically-based coupling relationships between climate systems. We 
then develop a method for creating paleoclimate networks using a range of proxy records. 
The resulting paleoclimate networks identify linkages at the global and regional scale, and 
demonstrate the transient nature of coupling relationships of the northwest Australian 
monsoon region throughout the Holocene. 


2.1 Data 


Table 1: Modern climate data 

Code 

Location 

Lat/Lon Bounds 

Season 

NWAusdjf 

Northwest Australia 

10-20°S; 115-140°E 

DJF 

NEAusdjf 

Northeast Australia 

10-20°S; 140-150°E 

DJF 

MCdjf 

Maritime Continent 

5°N-10°S; 90-150°E 

DJF 

IOdjf 

Western Indian Ocean 

0-15°S; 45-60°E 

DJF 

ISMjja 

Indian summer monsoon region 

5-25°N; 70-100°E 

JJA 

EASMj JA 

East Asian summer monsoon region 

10-20°N; 100-120°E 

JJA 

EEPdjf 

East Equatorial Pacific 

0-10°N; 230-250°E 

DJF 


Our main aim is to capture coupling relationships of the Holocene Australian summer 
monsoon, but we first test the suitability of complex networks to identify dynamically-based 
coupling relationships using modern climate data. Seasonal convective indices are con¬ 
structed using monthly values for 1948-2013 of mid-tropospheric (500mb) vertical velocity 
(w), a surrogate for convection (NCEP Reanalysis data provided by the NOAA/OAR/ESRL 
PSD, Boulder, Colorado, USA, from their web site at http://www.esrl.noaa.gov/psd/; 


Kalnay et al. 1996). In order to capture only coupling between deep convection, such as 


that associated with the monsoon circulation, we extract only three months of data from 
each year: December to February (DJF) or June to August (JJA), setting the values for 
the other nine months to zero (Table [I]). This data is averaged over the regions covering 
northwest Australia (NWAusdjf)j northeast Australia (NEAusdjf); the Maritime Con¬ 
tinent (MCdjf) ; the western Indian Ocean (IOdjf), the Indian summer monsoon region 
(ISMjja), the East Asian summer monsoon region (EASMjja), and the Eastern Equatorial 
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Figure 1: a) DJF 1981-2010 500mb u (NCEP Reanalysis, Kalnay et al., 1996). Also shown 
are the location of the boxes over which 500mb omega is averaged to produce modern 
convective indices (see Table [I]), and the location of proxies; b) As above, for JJA. 
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Table 2: Proxy records used in this analysis 


Code Location 

Lat/Lon 

Proxy Type 

Reference 

Average Time 





Step (years) 


F07 

Qunf Cave, Oman 

17.17°N, 54.30°E 

Speleothem <5 1S O 

Fleitmann et al., 2007 

7.7 

M14 

Lonar Lake, India 

19.98°N, 7G.51°E 

Multi-proxy 

Menzel et al., 2014 

18.8 

H08 

Heshang Cave, China 

30.45°N, 100.42°E 

Speleothem <5 18 0 

Hu et al., 2008 

7.8 

DOS 

Dongge Cave, China 

25.28°N, 108.08°E 

Speleothem <5 18 0 

Dykoski et al., 2005 

14.7 

Y07 

Lake Huguang Maar, China 

21.15°N, 110.28°E 

Ti concentration 

of lake sediment 

Yancheva et al., 2007 

0.8 

D10 

Sanbao Cave, China 

31.67°N, 110.43°E 

Speleothem <5 18 0 

Dong et al., 2010 

10.2 

G09 

Liang Luar Cave, Indonesia 

8.52°S, 120.43°E 

Speleothem <5 18 0 

Griffiths et al., 2009 

10.1 

D13 

Cave KNI-51, Australia 

15.30°S, 128.62°E 

Speleothem <5 18 0 

Denniston et al., 2013 

6.2 

M02 

Laguna Pallacocha, Ecuador 

2.77°S, 79.23°W 

Red colour intensity 
of lake sediment 

Moy et al., 2002 

0.8 

vB08 

Cueva del Tigre Perdido, Peru 

5.94°S, 77.31°W 

Speleothem <5 18 0 

van Breuklen et al., 2008 

19.4 

H01 

Cariaco Basin 

10.70°N, 65.17°W 

Ti concentration 
of marine sediment 

Haug et al., 2001 

5.6 


Pacific (EEPdjf)- Note that the use of convective indices prevents the incorporation of 
the East Asian winter monsoon in our analysis. The East Asian winter monsoon is char¬ 
acterised by northerly winds driven by the Siberian High, causing cold surges outflowing 
over the South China Sea. There is some related convective activity in southern China, 
but insufficient to be captured by a convective-based index. 

Following this, paleoclimate networks are produced for rolling 3,000 year windows at 
millennial intervals over the period 9,000 years BP to Present. We select proxy records 
(Table [2j) within the broad Indian Ocean-Pacific region according to high temporal reso¬ 
lution and low age uncertainty, as per Rehfeld and Kurths (2014). Although one prefers a 
database comprised of a single proxy for reasons of comparability, one is often constrained 


et al. 2007J Hu et al. 2008| Dykoski et al., 

2005f 

Dong et al. 

2010J Griffiths et al. 

2009| 

Denniston et al. 

2013b; 

van Breukelen et al. 

2008 

), titanium 

Yancheva et al., 2007 

Haug 


et al., 2001), sediment (Moy et al., 2002) and multi-proxy (Menzel et al., 2014) data sets. 


The IASM region is represented in the proxy record database by two speleothem records, 
G09 (Liang Luar, Flores; Griffiths et al. 2009) and D13 (Cave KNI-51, northwest Australia; 


Denniston et al., 2013b), both of which are interpreted as capturing monsoon precipitation 


trends and variation. The Chinese speleothem <5 18 0 records (D05, H08, D10) have each 
been interpreted as a proxy for precipitation changes driven by the East Asian summer 
monsoon, while the Lake Huguang Maar record (Yancheva et al. 2007) has been discussed 
in the context of the East Asian winter monsoon and coupled to the IASM region in the 
modern climate through cold surges. We also include two widely used proxy records: the 


titanium concentration series from the Cariaco basin (H01; Haug et al. 2001) has been 


cited in studies in the context of Holocene ITCZ positioning, and the Laguna Pallacocha 
sediment record from Peru (M02; Moy et al. 2002) is a very widely used proxy for changes 
in El Nino intensity and frequency over the last 12,000 years. 
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2.2 Constructing complex networks 


Estimating correlations between paleoclimate records is fraught with difficulty, and there¬ 
fore an intuitive qualitative curve-fitting approach is typically employed. We apply methods 
widely accepted by statistical physicists which have been successfully applied in the con- 


text of financial markets ( 

Zhuang et al. 

2014 

), solar activity ( 

Zou et al. 

2014 

), disease 

dynamics ( 

Zhang et al. 

2010 

Wu et ah, 

2015, 

Li et al., 

2015), and pigeon interactions in 

flight (Dieck Kattas et al., 

2012, 

Xu et al. 

2012 

). In a climate or paleoclimate context, 


one may envisage such a network as a number of nodes, each corresponding to the site of a 
climate or paleoclimate data set. If a statistically significant ‘similarity’ between two data 
sets is found, then an edge is drawn between the two nodes. More formally, for a database 
of n time series, denoted X % , we may describe the set of nodes as V = {vi : i E [n]}, and the 
set of edges is given by E = {ejj} where = 1 is Xj and Xj are found to be statistically 
significantly ‘similar’, and = 0 otherwise.We define similarity between two time series, 
Xi and Xj, by mutual information, a nonlinear, symmetric (and thus non-directional) mea¬ 
sure of how much information is shared between the two time series. Mutual information, 
I(Xi,Xj) is given by: 


I(X u Xj) = 


£ £ 

XiEXi Xj^Xj 


p(xi,Xj) log 


p{xj,Xj ) 
p(Xi)p(Xj] 


where p(xi) is the probability mass function of random variable Xi, and p(xi,Xj ) is the 
joint probability mass function of Xi and Xj. Note that, if Xi and Xj are independent, 
p(xi,Xj) = p(xi)p(xj), and hence mutual information is zero. If they are not independent, 
then the amount to which p{xi,Xj) differs from the product p{xi)p{xj) provides a mea¬ 
sure of the similarity of the two time series. We interpret this as a measure of coupling 
strength, with the information transfer between climate indices occurring through physi¬ 
cal atmospheric flows and pressure-driven teleconnections. We choose mutual information 
over linear cross-correlation due to the nonlinear nature of the relationship between pairs of 
proxy records, visible in scatterplots. Cross-correlation can produce spurious results in this 
situation (Kantz and Schreiber, 2003). We note the bias inherent in the Gaussian mutual 
information estimate, as demonstrated by Rehfeld and Kurths (2014). Irregular downsam¬ 
pling causes mutual information to be underestimated, and this bias increases rapidly when 
there are below 80-100 data points. Our analysis accounts for this: all the proxy records 


have more than 95 data points in any analysis window, except M14 (Lonar Lake, Menzel 


et ah, 2014) which has only 79 during the period 5,000-2,000 years BP and 86 during in 


the window 9,000-6,000 years BP. Rehfeld and Kurths (2014) suggest a bias-correction 
method for comparison with other similarity measures such as cross-correlation, but this 
is unnecessary here. Another similarity measure, the event synchronisation function, has 
been suggested to measure coupling between extreme events (Rehfeld and Kurths, 2014). 
However, this requires the use of only the data beyond, say, the 90% percentile, which 


7 

























































for the modern convective indices would provide only around 20 data points in each time 
series. As a result, it is unsuitable for this analysis. Prior to estimation, the raw time series 
data are detrended using a Gaussian high-pass filter with a bandwidth equal to half the 
analysis window. For the paleoclimate time series, this is equal to 1,500 years, and means 
we remove frequencies slower than 1 oscillation every 9.4kyrs. Only the non-zero data 
points in the modern convective indices are detrended and used in the subsequent analysis. 
Since the teleconnections between the regional monsoon systems often involve a time delay 
of up to 12 months, in the modern data analysis, we estimate mutual information over a 
-12 to +12 month window, and take the maximum value. 

Paleoclimate time series are often distributed along irregular time intervals due to 
sampling constraints. To account for this, a Gaussian kernel is used to ‘match’ data 
in paired paleoclimate time series. 


Rehfeld et al. 


(2011) demonstrate that this reduces 
bias in the resulting mutual information estimate compared to linear interpolation. We 
use the Matlab toolbox of Rehfeld and Kurths (2014) to produce estimates of Gaussian 
kernel weighted mutual information, Iq^Xi, Xj). This method does not produce symmetric 
estimates of Iq, but these asymmetric estimates do not imply directionality in the network, 
and are simply due to the unequal sampling rates of the two paleoclimate time series 
(Rehfeld et al., 2011). We therefore define: 


lG{Xi,Xj) = maxilciXi, Xj), I G {Xj, W)) 

We use a Monte Carlo approach to define statistically significant coupling relationships. 
For each modern or paleoclimate data set we generate a synthetic time series uncoupled 


to the others. Following Rehfeld et al. (2013), we use an autoregressive model with one 


lag, Brownian motion with drift, to model the modern data sets and all but one of the 
paleoclimate data sets. The parameters - linear drift and constant diffusion - are estimated 
from the observed time series through linear regression. This time series is initially regularly 
spaced, and we downsample according to the time steps of the original, observed data 
set. The Laguna Pallacocha record from Ecuador (M02, Moy et al. 2002) is not well 


suited to be modelled by Brownian motion. This time series is comprised of a number of 
large events which are registered well above a baseline level of near zero. We therefore 
introduce a Poisson process, to model the event time series defined by the 90% quantile in 
the Laguna Pallacocha record. This event time series is well approximated by a Poisson 
process (y 2 = 1.85, p = 10.12, at a 95% significance level). 

The synthetic time series are used to determine statistically significant coupling rela¬ 
tionships. Using the random time series models above, we generate synthetic time series 
corresponding to the observed data sets and create pairwise mutual information estimates. 
Repeating this 2000 times, we create 2000 randomly generated mutual information esti¬ 
mates for each pair of records. If the mutual information estimate calculated from the 
observed (modern or paleoclimate) data sets exceeds the 95th percentile of the randomly 
generated estimates, we claim there is a statistically significant coupling relationships be- 














tween the two records. Only connections which are identified as statistically significant are 
displayed in the network. 

Having constructed networks for the paleoclimate database (Table [ 2 ]) at 3,000 year 
windows throughout the last 9,000 years, we seek to evaluate changes in network density 
and structure. This may be attempted through a number of measures provided by graph 
theory (Newman, 2010). The degree, di, of a node, v^, describes the number of edges 
incident to the node, providing a description of how coupled the time series at Vi is to 
other records in the network. Similarly, the network average degree, d n , is given by: 


dn 


1 

n 


i =1 


This quantifies the total amount of coupling within the network. In addition to total net¬ 
work connectivity, we consider the degree distribution, the probability distribution of di 
across the network. This allows us to determine whether the modern climate or paleocli¬ 
mate records are all coupled to a similar degree, or whether a few records in particular are 
more dominant, potentially driving the broader monsoon network. In addition, we compare 
the observed degree distribution with the one from a random network where any two nodes 
are connected with probability p = |£ , |/|^t i |, where Er gives the total number of possible 
edges. Note that this produces a binomial degree distribution, taking as its parameters the 
number of nodes in the network and the number of connected edges as a fraction of the 
total number of possible edges. We then compare our observed degree distributions to this 
to identify any significant skew in the network connections. 


3 Testing complex networks using modern climate data 


Because the coupling relationships between modern regional climate systems are relatively 
well known, they provide a suitable control against which we evaluate the use of complex 
networks. Figure [2^, displays the known dynamical mechanisms by which the Indian, 
East Asian and Indonesian-Australian monsoons interact. Note that, with the Indian and 
East Asian summer monsoons active in boreal summer (JJA), and the East Asian winter 
monsoon and Indonesia-Australian summer monsoon occurring in austral summer (DJF), 
some of these interactions occur with a seasonal lag. 

In providing an explanation for the coupling relationships recognised in the ‘modern’ 
data, we initially appeal to Chiang’s (2009) framework for understanding the climate of 
the tropics. Chiang (2009) outlines two models of tropical circulation. Providing a basic 


mechanism, the climate of the topics is explained by the Hadley circulation and its response 
to seasonal heating. Here, the migration of the ITCZ is primarily responsible for the 
distribution and timing of precipitation across the tropics. However, he notes that this 
simple explanation is only sufficient in an aquaplanet setting, and that the existence and 
location of land masses introduces regionally. It is within this second conceptual model 
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Figure 2: a) Recognised interactions within the Asian-Australasian monsoon systems. We 
note the northward (JJA) and southward (DJF) positioning of the ITCZ, as well as i) the 
convective centre of the Bay of Bengal associated with the Indian summer monsoon, ii) 
northerly cold surges associated with the East Asian winter monsoon, iii) the Indo-Pacific 
Warm Pool, iv) the Nino 3.4 region, and v) the Madden-Julian Oscillation. For reference, 
the box denotes the northwest Australian monsoon region as defined for this study, b) 
Modern climate network. Nodes are located in the centre of the zonal averaging region 
given in Table [TJ Linked nodes are considered to be coupled at the 95% significance level. 












that the regional monsoons are explained, asymmetries recognised, and the role of ENSO 
incorporated. 

Figure[2jr depicts the ITCZ positioning in boreal and austral summer, with convergence- 
driven convective activity located over India, southern China and northern South America 
in JJA, while in DJF the ITCZ sits south of the equator, bringing convective activity to 
the Indonesian-Australian region and the Indian Ocean stretching from Indonesia to the 
north tip of Madagascar. Interactions between different regional features are also depicted. 
Outflows from the Bay of Bengal, associated with the Indian summer monsoon, are a key 
moisture source in the East Asian summer monsoon (Yihui and Chan, 2005). The Bay of 
Bengal is also the point of origin of convective centres which are displaced southwards over 


a number of months towards the Indonesian-Australian monsoon region (Meehl, 1987) 


Hung et al. (2004) demonstrate a correlation between the Indian and Australian monsoon 


regions, but describe a “communication gap” between the two systems, with heavy (weak) 
Indian summer monsoon precipitation followed by heavy (weak) Australian monsoon pre¬ 
cipitation. It is worth noting that this analysis combined the northwest and northeast 
of Australia into a single region, creating correlations with the El Nino-Southern Oscilla¬ 
tion (ENSO) which are likely due to the impact of ENSO on tropical northeast Australia. 
The relationship between ENSO and the Indian (Krishnamurthy and Goswami, 2000) and 
East Asian (Wang et al. 2000) summer monsoons is well established. Shifts in both the 
Walker and Hadley cells caused by warm (cool) sea surface temperature anomalies in the 
East Equatorial Pacific act to dampen (strengthen) the Indian summer monsoon, while 
El Nino events establish Rossby waves travelling towards China, setting up a region of 
anticyclonic circulation over the Philippine Sea and suppressing East Asian summer mon¬ 
soon convection. Finally, the East Asian winter monsoon establishes northerly winds over 
China, producing irregular low-level surges of cool air which travel southwards across the 
South China Sea and into the Indonesian sector (Chang et al., 1979). These surges are able 
to enhance convective activity, uplifting the already warm, moist air. A second outflow of 
cold air has been proposed to travel westward, flowing to the north of the Tibetan Plateau 
before being deflected southwards and across the Arabian Sea. Wang et al. (2012a) argue 
that this influx of cool air excites the MJO, thus acting as a secondary forcing mechanism 
on the Indonesian-Australian monsoon regime. 

The mutual information-based climate network (Figure [2^) captures many of the in¬ 
teractions outlined above. Due to the fact that the convective indices incorporate only 
boreal summer data from the East Asian region, we cannot capture any of the interactions 
involving the East Asian winter monsoon region. As noted in the methods section, this is 
intentional, as our cu-based indices cannot capture East Asian winter monsoon strength. 
The climate network, however, does miss two accepted coupling relationships. We expect 
a link to be observed between the east equatorial Pacific (EEPdjf) and northeast Aus¬ 
tralia (NEAusdjf)- There are two possible reasons for not capturing this in our analysis: 
either the link is missing because the signal cannot be distinguished through the noise of 
the climate time series, or the box size of the northeast Australian region (cf Figure [lj 
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might be unable to accurately capture regional precipitation using the NCEP/NCAR re¬ 
analysis. This is an issue which translates into the paleoclimate context directly, as many 
proxy records will not provide a ‘pure’ signal of monsoon-related precipitation, but rather 
capture a number of other climate and environmental changes. 

The second link missing in the modern climate network is between the Maritime 
Continent (MCdjf) and northwest Australia (NWAusdjf)- Given that the Indonesian- 
Australian summer monsoon extends across both regions, a coupling relationship between 


the two time series would be expected. However, Haylock and McBride (2001) examine 


summer rainfall measured at 63 stations across Indonesia, demonstrating limited spatial 
coherence across the region, with no single forcing mechanism or predictor of wet season 
precipitation. As such, averaging DJF convective activity across the region is unlikely to 
produce an index which can be interpreted easily in the context of coupling relationships 
with other convective indices. 

Despite the two missing network links, we have confidence in the ability of complex 
networks to capture coupling relationships between climate signals. The fact that there are 
no coupling relationships identified which do not correspond to any understood dynamical 
mechanisms further supports this, and validates the decision to produce networks based on 
mutual information. The spurious values which can arise when estimating cross-correlation 
between non-linear time series could lead to coupling relationships being identified which 
have no physical basis. 


4 Coupling relationships of the Australian summer monsoon 
over the last 9,000 years 

The Holocene combines, among other things, a period of changing solar insolation forcing, 
sea surface temperature feedbacks, and changes in land extent following the last deglacia¬ 
tion (cf. |Lambeck and Nakada 1990 Liu et al. 


2003, Jansen et al. 2008). We would 


therefore expect that complex networks may provide a formal demonstration of changing 
teleconnections throughout the Holocene. We demonstrated in section 3 that complex net¬ 
works can capture coupling relationships between modern climate signals. We now present 
paleoclimate networks at 3,000 year windows, every millennium, from 9,000 yrs BP to 
Present. The network is populated by the database of proxy records from the Asian- 
Pacific region (Table [2]). Each edge in the network (Figure [3]) identifies a statistically 
significant amount of information shared between two paleoclimate proxy records. We in¬ 
terpret these edges as dynamically-based coupling relationships between regional climate 
systems. With network edges identifying significant coupling relationships, the total degree 
of connectivity within the network should provide insight into the collective behaviour of 
regional monsoon systems during the Holocene. However, as our focus is on the Australian 
summer monsoon, we emphasise nodes G09 and D13 (corresponding to Liang Luar, Flo¬ 
res, Indonesia; Griffiths et al., 2009, and Cave KNI-51, northwest Australia; Denniston 



































|Figure 3: Paleoclimate networks for a: 9,1^)0-6,000 yrs BP, b: 8,000-5,000 yrs BP, ~cT] 

7,000-4,000 yrs BP, d: 6,000-3,000 yrs BP. Linked nodes are coupled at the 95% level. 






















Figure 3 (Cont.): e: 5,000-2,000 yrs BP, /: 4,000-1,000 yrs BP, g : 3,000-0 yrs BP. 
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et al. 2013b) to enable us to draw inferences regarding the coupling relationships of the 


Indonesian-Australian monsoon regime. 


4.1 Overall trends in network relationships 

The paleoclimate networks are observed to grow increasingly connected from 9,000-6,000 
yrs BP to 6,000-3,000 yrs BP, and then decline steadily into the latest Holocene (Fig¬ 
ure J3]). In particular, the nodes corresponding to the records H01 (Cariaco Basin; Haug 


et al. 2001), M02 (Laguna Pallacocha, Ecuador; Moy et al. 2002| ), vB08 (Cueva del Tigre 
Perdido, Peru; van Breukelen et al., 2008), and Y07 (Lake Huguang Maar, China; Yancheva 


et al. 


2007) display high degree throughout the mid-Holocene, while only M02 is seen to 
be highly coupled in the later Holocene. We quantify total network connectivity using 
the average network degree. This highlights a trend of increasing degree until 6,000-3,000 
yrs BP, followed by a return to nodes of lower degree (Figure |4^i). To determine if the 
coupling relationships are evenly distributed across the network we compare the observed 
degree distribution for each 3,000 year window with the one of a random graph. This 
random graph has the same number of nodes and edges as the observed network and a 
binomial degree distribution. The small dataset accounts for the large confidence error 
bands, which mean that we cannot state much, with confidence, about the difference be¬ 
tween our observed networks and a random graph. However, during the mid-Holocene the 
degree distributions resemble that of a random graph which sits within the 90% confidence 
bands (Figures [4j: - e). This means that although the network has a higher degree during 
this period, coupling relationships are spread somewhat evenly across the region - rather 
a situation where one proxy record dominates the network, instead there is a relatively 
equal amount of coherence between the regional monsoon systems. In the late Holocene 
(Figures [4^ - h) there is increasing skewness, with a statistically significant deviation for 
degree greater than 7 from the distribution from a random graph. This demonstrates 
that although most proxy records have few coupling relationships with other records, a 
small number of records ‘dominate’ the network, playing a more critical role. The Laguna 
Pallacocha (M02) record stands out here as being very highly connected, with coupling 
relationships identified with nearly every other record in the network. 


4.2 Coupling relationships of the I ASM region 

The proxy records within the Indonesian-Australian monsoon region are G09 and D13, 


relationship between the two records is identified in the paleoclimate networks (Figure [3]). 
Possible explanations for this include the rising sea levels experienced during the deglacia¬ 
tion and into the mid-Holocene, and by the modern-day climatology which displays large 


speleothem d 18 0 records located in Indonesia and northwest Australia respectively (Grif- 


fiths et al. 

2009 

Denniston et al. 

2013b) 


ing local precipitation, and tied directly to the IASM circulation. Despite this, no coupling 
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Figure 4: a) Network average degree; b-h) Degree distributions for the observed paleocli- 
mate network (black), with 90% upper and lower confidence intervals (dotted lines), and 
that of a random graph with the same number of edges (grey). 
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spatial variations in monsoon precipitation across the region. To investigate this, we draw 
upon other proxy records from the region, unsuitable for our network analysis, but able to 
assist in a qualitative interpretation. 

With the progressive increase in sea level from c. -130 m over the 20 000 years (sum¬ 
mary in Murray-Wallace and Woodroffe, 2014), land extent over the IASM region was 


significantly modified (Voris 2000). During our period of interest, records suggest a sea 


level of around 20m below PMSL at around 9,000 years BP, increasing to present day levels 
by about 7,000 years BP ( Lewis et al.[ |2013{ |Murray- Wallace and Woodroffe 2014). Pa- 
leoclimate proxy studies from the IASM region suggest a broad-scale response to sea level 
changes in the early to mid-Holocene (Figures^, b, c). Griffiths et al. (2009) interpret the 


strengthening precipitation signal from around 11,000 to 6,000 years BP in the Liang Luar 
(G10) speleothem record as evidence of a monsoon response to the flooding of the Sunda 
shelf. This is also noted by Denniston et al. (2013b) in the Cave KNI-51 record (D13), 
with strengthening from 9,000 to 7,000 years BP possibly tied to Sahul and Sunda shelf 
flooding. A similar strengthening of the monsoon is also recognised in other speleothem 
records: Gunung Buda, north Borneo (Figure [5fc,|Partin et al. 2007) and Ball Gown Cave, 


northwest Australia ( |Denniston et ah] |2013a| ) . These trends represent a multi-millennial 
scale response to sea level change, with monsoon strengthening observed over some 5,000 


years. 


The Indonesian (G09) and northwest Australian (D13) records remain uncoupled in 
the later networks, which clearly cannot be attributed to sea level rise. Instead, we look 


to the spatial variations observed in the modern IASM circulation (Haylock and McBride 
2001). Heterogeneity in monsoon-related precipitation across the IASM region has been 


attributed to possible local relief or island controls (e.g. Moron et al. 2009) as well as 
varying interactions with the Indo-Pacific Warm Pool and sensitivity to interannual ITCZ 
positioning (e.g. Partin et ah , 2007). Proxy records available in the late Holocene - Flores 


Indonesia (G09, Griffiths et ah, 2009), Cave KNI-51, northwest Australia (D13, Denniston 


et ah, 2013b|), Gunung Buda, northern Borneo ([Partin et ah 2007), and Lombok Basin 


near Sumba, Indonesia ( |Steinke et al.[ |2014| ) - indicate no coherence in monsoon strength 
at millennial to centennial time scales (Figure [5]). |Denniston et al. (2013b) consider their 
Cave KNI-51 record to be anti-phased with the Flores record in the late Holocene, and 
note that the weakening of the monsoon observed in the northwest Australian record from 
4,000 to 1,500 years BP is much less visible in the Indonesian record. The Lombok Basin 


record (Steinke et ah 2014) has been interpreted as displaying a very different chronology, 


with an abrupt increase in precipitation at around 2,800 years BP, with this stronger 
monsoon regime continuing until 1,700 years BP. An explanation for the difference in sub- 
millennial scale events in the northwest Australian and Indonesian records during the late 
Holocene was considered by Steinke et ah (2014). They suggested that modelling results 
of the solar minimum at around 2,800 years BP demonstrated varied responses across 
the IASM region, with strengthening across Indonesia, weakening in northern Australia, 
and a somewhat neutral response in Borneo. It is therefore reasonable that short-term 
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Figure 5: Holocene proxy records from the Indonesian-Australian monsoon region. 


speleothem <5 18 0 record from Cave KNI-51, northwest Australia (D13; Denniston et al. 


c) speleothem d 18 0 record from Gunung Buda, Borneo (Partin et al., 2007), d) sediment 


record from offshore Surnba Island, Indonesia (Steinke et al., 2014). The arrows denote in 


2013b 

), b) speleothem d 18 0 record from Flores, Indonesia (G09; 

Griffiths et al., 

2009) 


the direction of enhanced monsoon strength in the proxy record interpretation. 


fluctuations in the IASM paleoclimate proxy records might be uncoupled due to spatially 
heterogeneous responses to external forcing across the region, just as can be observed in 
the present day climatology. 

Despite the differences between the Liang Luar (G09; Griffiths et al., 2009) and Cave 
KNI-51 (D13; Denniston et al.[ |2013b| ) records, they display similar coupling relationships 
with other, non-IASM proxies (Figure[3]). From 9,000-6,000 yrs BP to 4,000-1,000 yrs BP, 
these two records are coupled (with some inconsistency over time) to number of records 
forming a belt from the Arabian Pennisula to South America: Qunf Cave, Oman (F07 
Fleitmann et al., 2007|), Lonar Lake, India (M14; Menzel et al., 2014|), Lake Huguang Maar 


southern China (Y07; |Yancheva et al. 2007|) , Cueva del Perdido, Peru (vB08; van Breukelen 


et al. 2008), and Cariaco Basin (H01; Haug et al. 2001). There is, additionally, some 


coupling between records in the Indonesian-Australian and East Asian summer monsoon 
regions during 6,000-3,000 yrs BP and 5,000-2,000 yrs BP. However, in the latest Holocene 
(3,000-0 yrs BP), the only record which both Liang Luar (G09; Griffiths et al. 2009) and 
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Cave KNI-51 (D13; Denniston et al. 

2013b 

) are connected to is the Laguna Pallacocha 

record from Ecuador (M02; |Moy et al. 

2002 

)• 


4.3 Interpretation of IASM coupling relationships of the last 9,000 years 


The paleoclimate proxy records from Oman, India, southern China, Peru and Cariaco Basin 
have each been interpreted in the context of ITCZ positioning. The titanium concentration 
time series from Cariaco basin (HOI; Haug et al., 2001| ) and the speleothem oxygen isotope 
composition in Cueva del Tigre Perdido, Peru (vB08 


van Breukelen et al. 2008) 


are 


dependent on the position of the ITCZ during boreal winter. The Lake Huguang Maar 
sediment record in southern China (Y07) is interpreted as a proxy for the East Asian 
winter monsoon strength, determined over millennial scales by the position of the ITCZ 
(Yancheva et al. 2007). The bioclastic climate index from Lonar Lake, central India (M14) 


is not explicitly interpreted as an indicator of ITCZ positioning, but the authors tie early 


Holocene Indian summer monsoon strength to ITCZ migration (Menzel et al. 2014 Prasad 


et al. 2014). Finally, Qunf cave in southern Oman (F07; Fleitmann et al., 2007) is situated 
at the southern edge of the modern position of the ITCZ in boreal summer. A northwards 
shift of the ITCZ causes southwesterly flow associated with the Indian summer monsoon to 
extend over the lower tip of the Arabian Peninsula, lifting the local temperature inversion 
and triggering deep convective precipitation. As the ITCZ shifts south, however, its role in 
connecting these two regions ceases. The coupling relationships identified between F09 and 
M14 from 9,000 to 4,000 years BP (Figure [3^,-c) therefore support a northerly positioning of 
the ITCZ, as expected under Northern Hemisphere precessional bias. At the same time, the 
monsoon region of Indonesia and northwest Australia experiences a weak but strengthening 
monsoon, and the density of coupling relationships with TTCZ-proxies’ and the East Asian 
summer monsoon proxies suggests an Indonesian-Australian monsoon modulated by global- 
scale forcing. 


The northwest Australian proxy from Cave KNI-51 (D13; Denniston et al., 2013b) 


displays the most connections with other proxy records during the period 6,000-3,000 yrs 
BP. This coincides with a period of dense coupling across the full network. Inspection 
of the raw data (Figure [5^,) shows this period to be transitional, with a step change in 
d 18 0 values at around 4,000 years BP. This would be expected with a shift to Southern 
Hemisphere precessional bias, and thus the paleoclimate network seems to capture this 
regional response to changes in Milankovich forcing. 

From 5,000-2,000 years BP onwards, the number of coupling relationships between 
the IASM proxies and the ITCZ-related and Chinese proxy records begins to decline. 
Widespread network connectivity also decreases, and the absence of coupling between F07 


(Qunf Cave, Oman; Fleitmann et al. 2007) and M14 (Lonar Lake, India; Menzel et al. 


2014| ) may indicate that the networks have captured the southward progression of the 
ITCZ. These two coincident trends suggest that global-scale synchronicity declines in the 
later Holocene because of decreasing direct insolation in the Northern Hemisphere. This 
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dampens the strength of the Indian and East Asian monsoon, thus weakening interhemi- 
spheric coupling relationships with the IASM region. As such, a model emerges whereby 
Milankovich insolation forcing acts as a control on regional monsoon strength not only 
through direct radiative forcing, but also indirectly, by modulating the strength of cou¬ 
pling relationships between regional monsoon systems. 

This model of indirect insolation control agrees with proxy records from the Indonesian- 
Australian monsoon region. Denniston et al. (2013b) note that although precession and 
tilt favour the Southern Hemisphere following 6,000 years BP, precipitation over northwest 
Australia is observed to decline. Liu et al. (2003) present sea surface temperature feedbacks 
as an explanatory mechanism for the counterintuitive response to precessional forcing over 
northwest Australia, while Wyrwoll et al. (2007) demonstrate that tilt as well as precession 
plays a critical role in determining monsoon precipitation over northwest Australia. Using 
a coupled ocean-atmospheric model they found that high tilt lead to enhanced monsoon 
precipitation, even under a Northern Hemisphere precession bias. In fact, the simulation 
results display changes in interhemispheric outflows to the Southern Hemisphere between 
different precession and tilt scenarios (Wyrwoll et al., 2007). The observed weakening in 
coupling relationships between the IASM region and the Northern Hemisphere into the later 
Holocene suggests that the transient nature of teleconnections between regional climate 
systems may have played a critical role in determining the response of the Indonesian- 
Australian monsoon to Milankovich insolation forcing. 

By the latest Holocene (3,000-0 years BP) the Laguna Pallacocha record in Peru (M02; 
Moy et al., 2002) is the only record coupled to both the IASM proxy time series. Be¬ 


cause this time series is not suitably represented by Brownian motion (unlike, say, the 
<5 ls O speleothem time series), we have defined statistical significance based on simula¬ 
tions of a surrogate Poisson process. This Poisson process is derived from an ‘event’ time 
series created from 90th percentile events within the M02 record, and is well approxi¬ 
mated by a Poisson process. Mutual information, as a nonlinear measure of similarity 
between time series, should be able to quantify the coupling relationships of the Laguna 
Pallacocha red intensity index (M02; |Moy et al., 2002). However, the mutual information 
estimates calculated from these surrogate Poisson time series will not correspond directly 
with estimates calculated using the original M02 record. Replicating this analysis using 
a bootstrapping approach to simulate the M02 record returned the same network struc¬ 
tures, indicating that this method is suitable. However, we still proceed with caution when 
drawing conclusions regarding coupling relationships with the Laguna Pallacocha record. 
In the latest Holocene, dense coupling with the M02 record is observed, including both the 
IASM records. Moy et al. (2002) interpret the events recorded in the red intensity index 
as indicative of moderate-to-strong El Nino events. They therefore argue that the record 
demonstrates an increase in these events from around 3,000 years BP until 1,200 years BP. 
Other marine and terrestrial proxy records from the Pacific region indicate similar trends - 
lake sediment records from Galapagos (Conroy et al., 2008), and proxy records of upwelling 
on the Peruvian margin (Rein et al. 2005) and in the Panama basin (Cabarcos et al. 2014). 
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However, at sub-millennial timescales, the timing of “spikes” in ENSO activity does not 
match between records, and there exists a further set of marine records which suggest 
slight different scenarios - A sediment core from the Peruvian margin suggests an increase 


in both El Nino and La Nina strength from 3,000 years BP onwards (Makou et ah, 2010) 


and Koutavas et al. (2006) argue that although the late Holocene experienced an increase 
in El Nino activity, El Nino strength was not abnormally high, but simply increasing after a 
period of La Nina-like conditions in the early to mid-Holocene. In the Australian summer 
monsoon context, a number of authors have suggested that enhanced ENSO activity in 
the latest Holocene may have contributed to a period of aridity evidenced in speleothem 


and paleoenvironmental records from northwest Australia (Shulmeister and Lees 1995 


McGowan et al. 

2012 

Denniston et al. 

2013b) 


However, it has been recognised for a long 
time that there is no apparent impact of ENSO on monsoon precipitation over northwest 
Australia (McBride and Nicholls, 1983). This is also evident in precipitation records, with 
the exception of the 1982- 1983 El Nino event, during which the Southern Oscillation Index 
reached a record low of -33 (Bureau of Meteorology, 2014a). During this event, there was a 
clear reduction in Australian summer monsoon rainfall. There is therefore the possibility of 
an ENSO-northwest Australian teleconnection, although only under an anomalously large 
shift in the Walker circulation. The paleoclimate networks provided here suggest, albeit 
tentatively, that a teleconnection between ENSO events and IASM precipitation may well 
have been present during the latest Holocene. However, given that the Laguna Pallacocha 
record is, itself, dependent on teleconnections with the central Nino 3.4 region, we make 
no firmer claims, as the relationship between the Peruvian and ENSO variability may also 
be transient during the later Holocene. 

In summary, by using the verified complex networks method to identify significant cou¬ 
pling relationships between paleoclimate proxy records, we are able to make the following 
comments: 


1. The observed coupling relationships between proxy records from the Indian, East 
Asian and Indonesian Australian monsoon regions highlight the validity of the ‘global 
monsoon’ model. 


2. However, these coupling relationships are observed to change over the last 9,000 years. 
In particular, a global-scale transition is observed, whereby the paleoclimate network 
decreases in connectivity towards the later Holocene. This may indicate reduced 
global-scale coupling caused by Southern Hemisphere precessional bias, causing a 
southward shift in the ITCZ. 


3. In the Australian summer monsoon context, connections are observed with a num¬ 
ber of Northern Hemispheric proxy records in the early to mid-Holocene. This 
corresponds with an understanding of interhemispheric flows modulating monsoon 
strength over northwest Australia (e.g. Liu et al., 2003). In the latest Holocene, 


these coupling relationships are lost and the Laguna Pallacocha record (M02; Moy 






























et al. 2002) is connected instead. This raises the possibility of an ENSO teleconnec¬ 


tion to northwest Australia during this period. This shows some alignment with the 
dual-model of tropical climate set out by Chiang (2009), with global scale controls 
(i.e. ITCZ) and those derived from regional differences (i.e. ENSO) both playing a 
role in determining regional monsoon strength throughout the last 9,000 years. 


5 Conclusions 

Identifying potential coupling relationships between climate systems using paleoclimate 
proxy records is typically a qualitative process. Here, we demonstrate the efficacy of com¬ 
plex networks to identify coupling relationships and teleconnections which correspond to 
known dynamical mechanisms. Extending this method to a multi-proxy database of pa¬ 
leoclimate time series, we are able to draw conclusions regarding the nature of coupling 
across the Asian-Australasian monsoon region. Our results recognise an element of the 
global monsoon concept, with regional monsoons displaying some degree of coupling over 
the period. The global paleomonsoon model, however, does not adequately represent the 
transient nature of coupling relationships, while our findings demonstrate a strengthening 
of coupling relationships across the broad Asian-Australasian monsoon regions during the 
mid-Holocene, followed by a tendency to reduced coupling in the later Holocene. Our 
findings at this stage are preliminary, and dependent on the availability of suitable proxy 
datasets, but we envisage that once more datasets become available, a stronger case can 
be made. In the context of the Australian summer monsoon, we observe coupling rela¬ 
tionships to other low latitude regions throughout the Holocene. While we offer tentative 
explanations for these, we stress that the observed links are unable to tell us about the 
underlying mechanism. We can, however, state with confidence that the networks demon¬ 
strate effectively that coupling relationships between the northwest Australian monsoon 
region and other regional climate systems were transient over the past 9,000 years. Given 
the coupling relationships observed, the next step is to ask why they exist. We recommend 
further research in the context of available model simulations, and in developing methods to 
assess the direction of information flow, in an effort to identify the underlying mechanisms 
determining and driving the coupling relationships. 
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